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Abstract 

A computer simulation method will be discussed that provides for equivalent 
simulation accuracy, but that exhibits significantly lower CPU running time per 
bias point compared to other teclmiques. This new method is applied to a taesh 
point field as is cus:oaary Jn numerical integration (NI) techniques. The 
assumption of a linear approximation for the dependent variable, which is 
typically used in the finite difference a^d finite element NI methods, is not 
required. Instead, the set of device transport equations is applied tOj and the 
closed-form solutions obtained for, each mesh point. The mesh point field is 
generated so that the coefficients in the set of transport equations exhibit 
small changes between adjacent mesh points. In contrast to the NI linear 
approximation, the closed-fora solutions more accurately represent the physical 
system and the device physics incorporated in the transport equations. 

Application of this method to high-efficiency silicon solar cells is 
described; and the met! >d by which Auger recombination, arcbipolar consider- 
ations, built-in and induced electric fields, bandgap narrowing carrier con- 
finement, and carrier dif fusivities are treated. Bandgap narrowing has been 
investigated using Fermi-Dirac statistics, and these results show that bandgap 
narrowing is more pronounced and that it is temperatu-e-dependent in contrast to 
the results based on Boltzmann statistics. It is al so suggested that carrier 
diffusivitv relationships that apply to degenerate materials in thermal equili- 
brium may also be applicable to regions in which high injection exists even in 
nondegenerate material. 

Imposing the appropriate boundary conditions on the u id-form solutions 
results in a ser of equations which require simultaneou so ..it ion. This results 
in obtaining the solution of all constants of integration, r rom which, in 
principle, all cell characteristics may be derived. It has been demonstrated 
that recursion relationship., exist between the constants of integration. Trial 
or “guess’ 1 solutions are not required in this new method for devices operating 
at any injection level, because the closed-form soli. -tons obtained at each mesh 
point, in fact, fulfull this cole . This carries over to those devices operating 
at high injection levels., but the inclusion of Auger recombination introduces 
nonlinear terms in the continuity equations, and special attention must be 
devoted to satisfying Poisson’s equation. Und r these conditions, an initial 
estimace must be made of the value of the independent variable for inclusion in 
the atinuity equations at th_ initial mesh point. An iterative procedure is 
then used to obtain a consistent solution. 
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INTRODUCTION 

Computer modeling simulations have been shown to be very useful in the 
development of semiconductor devices in those cases where the simulation is an 
accurate representation of the physical device. However, to be rn effective aid 
to the experimentalist and to become an equal partner in the technologies used 
in device development, it may be required to operate the computer program 
frequently each day in an active developmental program. For frequent use, as 
required in solar cell development, CPU costs must be low. Moreover, low CPU 
cost allows for engaging in computer experiments, which can be made to be a very 
useful and powerful technique. 

Computer modeling using numerical integration (NI) methods in Si device 
technology have usually shown fair-to-good agreement with experimental data. 
However, CPU costs for the execution of computer programs that are based on 
numerical integration methods are prohibitively high for their use as a labora- 
tory or manufacturing tool [!}. The number cf bias points that are required to 
study optimized device designs usually exceeds 5,000 runs. Similarly, a com- 
prehensive study involving device structures or new t-prs of devices may exceed 
10,000 runs. Increased CPU cost results if convergent problems arise. In most 
cases, the cost of such studies, for the benefits gained, may not be attractive. 

Simulation accuracy is determined by both the accuracy of the algorithm/ 
analytical method representing the device transport equations, and the accuracy 
of the phenomena submodels in representing the corresponding experimental data 
related to material properties. For most efficient use of the CPU, the accuracy 
of the algorithm/analytical method and of the phenomena submodels should be 
commensurate. For example, even if the algorithm/analytical method accurately 
represents the device, simulation results say not agree with experimental data 
if the phenomena submodels are accurately represente The reverse is also 
true. In solar cells, the phenomena submodels that produce first-order effects 
in terminal cnaracteristics are: abs' *ption curve, built-in and induced electric 

fields, bandgap, lifetime, mobilities, dif fusivities , photoexcited carrier 
concentrations, surface recombination velocities, junction transport, etc. The 
representation of the phenomena submodels must take on an importance equal to 
the analytical method used to represent the system. 

Under some operating conditions and for a number of solar cell structures, 
two-dimensional modeling may be required to obtain improved agreement between 
simulation results and experimental data. Although the results presented above 
apply generally to one- and two-dimensional modeling, the CPU execution time is 
signii icantly greater for two-dimensional simulations. 

In this paper, a new method is described which has been used to simulate 
semiconductor device characteristics. Although this method share, similarities 
with some aspects of HI methods, it differs markedly in other aspects. 

Abbreviated forms of analytical relationships representing the solution of solar 
cell transport equations that are obtained using this new method are presented 
and discussed. The similarities and differences between the methods are also 
discussed. 
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Simulation methods 


In this section, the transport equations are discussed as they relate to 
the simplifications and approximations that are required to be made in order to 
use HI meLhods and the new method, which is the application of closed-form 
solutions applied to a mesh p^int field. The finite difference (FD) and finite 
element (FE) methods, which are considered NI methods, are briefly discussed 
because they are most similar to this new method. An outline of these methods 
is presented and the procedures are developed to apply them. 


Numerical Integration Methods 

The most commonly used transport model in semiconductor devices is due to 
Van Roosbroeck and is represented by the set of equations: 


|f * | v * J n + (G-R) 
o t q n 

(1) 

If = - ~ V • J + (G-R) 
ot q p 

(2) 

J = ay nE + qD Vn , 

n ‘ n n 

(3) 

J = qy pi - qD Vp , 

P P P 

(4) 

V • E = 3 (p-n + N -N ) 
e DA 

(5) 


The above set of relationships is applicable to nondegenerate and degenerate 
materials, low and high injection levels, dc and ac operation, and to meat 
semiconductor device structures. 

The Van Roosbroeck transport equations comprise a system of coupled partial 
differential equations that describe a semiconductor carrier concentration and 
current density in position and time. The net balance of generation sources and 
recombination sinks of electrons and holes are described by the respective 
continuity equations, and the Poisson equation describes the electric field 
distribution that is produced by the charge distributior within the semiconductor. 
For solar cells, the steady-state condition is assumed, which greatly simplifies 
the transport set because dn/3t and 3p/H vanish. Device phenomena submodels 
may be added to tne transport equations to accurately describe a variety of 
carrier dynamical and other internal physical processes. 

However, the Van Roosbroeck equations must be signif icantl> modified i s . 
order to describe the effects arising from velocity overshoot, ballistic trans- 
port, and very thin surface layers [3. . In solar cells, only th- latter may 
need consideration. The phenomenon ir thin n-type surface layers, where an 
oxide charge insulator (OCI) produces an electron accumulation, may introduce 
sidebands within the conduction ar.J where each sideband corresponds to a 
quantized level for electron transport in a direction normal to the surface. 
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Solution of the transport equations is more dif f icult for higher order 
models; i.^., time-dependent requirements, and two- and three-dimensional 
geometries [2-11]. Moreover, this difficulty is brought into sharp focus in 
two-dimensional modeling when comprehensive phenomena submodels are considered 
[ 2 , 12 ]. 

In comprenensive representations, the complexity of the Van Roosbroeck 
transport equations does not allow for closed-form solutions which accurately 
simulate semiconductor devices. Accurate solution of the transport equations 
requires a method that simplifies these equations. The most popular approach is 
to divide the device structure into small parts [2-15]. In one-, two-, and 
three-dimensional solutions, these small parts are defined as infinite slabs, 
areas, and volumes, respectively. Each of these small parts is assigned a mesh 
point, which is identified by one, two, or three indices corresponding to a one-, 
two-, or three-dimensional representation of the device. These small parts must 
be made sufficiently small so that ail dependent variables in the transport 
equations exhibit small changes In value between adjacent mesh points [2, 4, 5, 8, 
15]. In addition, the coefficients in the transport equations also exhibit 
small or negligible changes between adjacent mesh points [15]. If either of 
these conditions is not met, a solution is not obtained because convergence does 
not occur. NI methods and the closed-form solutions may be applied to mesh 
point fields. 


In applying the FD method, all derivatives are replaced by finite differences 
be veen discrete points in an active domain in the interior of the structure. 

The residue of the newly established difference equation is set to zero at each 
me; n point. Thus, the differential equations are transformed into difference 
op rator equations. The value of the dependent variable is determined at each 
m sh ooint from the set of equations obtained [2,4]. 

For examp' e, the differential equation to determine the temperature distri- 
bution is given by [4] 





0 


( 6 ) 


The difference operation is represented by 



i+2 


- T 


i+1 


Ax 



(7) 


where T^, T^ + ^, and are the values of the dependent variable, T, at the 

mesh points i, i+1, and i+2, and Ax is the independent variable representing the 
separation of the mesh points i and i+3 , and i+1 and i+2, assumed to be equal. 

A cor * 'ponding equation is constructed at each mesh point. The resultant set 
o f equations requires the simultaneous solution for T^, T„, — , , — . Thus, 

v-1 >e of the dependent variable, T, is obtained o nly the mesh points 1, 

a. --. The value for T in the region between the mesh points is not 
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obtained explicitly and may only be approximated by one of the established 
methods. Typically, a linear approximation is used. 

In the FE method, a typical assumption made is that the dependent variable 
is a linear function, where, for the system discussed above involving temperature 

[4], 

^i+1 ~ ^i 

T * T _r + (x - x ) . (8) 

* i+1 i 1 

This relationship is substituted into an equation, called the "functional," that 
describes the system to obtain a set of equations. Simultaneous solution of 
this set results in discrete values of the dependent variable at each of the 
mesh points. Between the mesh points, the dependent variable is gr^erned by a 
linear relationship similar to that represented in Eq. (8). 

The full numerical solution of partial dif ferentic. equations as described 
above, which is applied to a physical system and which describes all regions of 
this system in a unified manner, was first suggested by Gummel [5] . In this 
work, Gunmel applied his method to a one-dimensional bipolar transistor. 
Subsequently, the method was further developed by DeMari [6,7] who applied it co 
p-n junctions, and by Scharfetter and Gummel [8] to IMP ATT diodes. Although 
this initial work was confined to one- dimensional structures, the Gummel approach 
has also been adapted to two-dimensional modeling [3,9-12]. However, the two- 
dimensional algorithm requires excessive CPU execution time [3,12], 

Phenomena submodels may be imposed on the coupled nonlinear partial 
differential equations which may result in nonlinear transport equations; i.e. , 

2 2 

band-to-band Auger recombination introduces a term involving n or p , and the 
electric field in the quasi-neutral region depends on the injection level. In 
NI numerical methods, the transport relationships are linearized. Two schemes 
have been used and both require initial "guess" solutions. The decoupled method 
proposed by G umm el [5] is to assume the coupling is weak so that the equations 
are solved serially. While the method is not difficult to implement, it fails 
to give accurate results for highly nonlinear systems. The coupled method, 
proposed by Hach*’ ,: ;l et ^1. [14], solves the transport equations simultaneously. 
However, this is mo»*e difficult to implement. The implementation of the 
algorithm becomes complex when recombination-generation and electric field drif 
terms are included [4]. In addition, this results in significantly higher CPU 
execution time and increased main memory requirements. 

As discussed above, .Jog NI methods results in values for the dependent 
variables, such as the c * _er concentrations and electric field, at discrete 
points in the semiconductor; i.e., at the mesh points. Thus, the -ontinuum of 
the dependent variables, described by the transport equations, is transformed 
into a discontinuous or discrete set cf values representing these variables at 
each of the mesh points. In contrast, the method that is proposed in this paper 
does not resort to this transformation. 
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Closed-ronn Solutions Applied to a Mesh Point Field 


The major point of departure between the algorithms of the new method and 
NI methods is that the new method uses a closed-form solution which is applied 
to a mesh point field that defines the system in space and in time [15]. This 
requires the application of the transport equations to each mesh point. In 
order to obtain a closed-form solution to the transport equations, the mesh 
point field is generated so that the equations accurately represent each mesh 
point. The closed-form relationships represent analytical solutions and result 
in a continuum of values for the dependent variables in the region between 
adjacent mesh points- Analytical solutions of adjacent mesh points also provide 
a continuum of values for the dependent variables in their corresponding regions. 
Analytical solutions of neighboring mesh points satisfy the usual boundary 
conditions that are demanded in semiconductor devices in accordance with the 
electronic/optical model under consideration. The analytical solutions repre- 
senting the dependent variables at each mesh point contain constants of integra- 
tion, which are determined by imposing the boundary conditions on each of the 
solutions. 


A procedure has been established in applying the closed-form solutions to a 
mesh point field. The procedure is as follows [15]: (1) establish the elec- 

tronic/optical model of the solar cell within the generic transport equations; 
(2) impose the phenomena submodels, represented in analytical or tabular form, 
on the transport equations; (3) develop mesh point field distribution of order f 
which reduces the complexity of, and makes the coefficients that are present in, 
the continuity equations constant or nearly constant, so that a closed-form 
solution is obtained at each mesh point with minimum restrictions; (4) establish 
2f-boundary conditions on the mesh point field; (5) apply the 2f-boundary 
conditions on the f-closed-form solutions; and (6) solve the resultant 2f~ 
equations simultaneously to obtain the 2f-constants of integration through 
recursion relationships. In principle, the transport equations are uniquely and 
completely solved after the relationships for the 2f-constants of integration 
are obtained. Electron current is calculated at the depletion region edge; the 
hole current is determined by a similar relationship. Adding the electron and 
hole currents results in the J-V relationship from which most of the terminal 
characteristics are recovered. 


The well-established solar cell electronic/optical model is imposed on the 
transport equations; i.e., optical pair generation of carriers, and minority 
carrier collection and transport across a p-n junction. This is followed by 
imposing the phenomena submodels which are subsequently discussed. A mesh point 
distribution is established so that the electric field, lifetime, mobility, 
diffusivity, absorption coefficient, and bandgap exhibit small changes between 
adjacent mesh points. Figure 1 show' the subdivision of a one-dimensional solar 
cell structure. The n-region (0, X^) and the p-region (X^, X,.) are subdivided 

into £ - and f -mesh points, respectively. To illustrate the method, the p- 

region transport solution is discussed below. The electron continuity equation 
is represented by 


So' t,N ° eKp (- Jo adx ') 


d^r. (x) 


d\ + D 




n 


dx 


+ U E 
n p 


an (x) n (*.) - n 
pe p __ 

dx i 


- 2 ° = 0 , 


(9) 
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Figure 1. Discretization of 1 -dimensional solar cell structure. 



where the Poisson equation takes its usual form. The symbols have their usual 

definitions, where n and n represent the total and photoexcited electron 

P pe 

concentrations, respectively. Implicit in Eq. (.9) is that the subdivision 
results in a condition in which the slope of the electric field, dE^/dx, is 

small and may be neglected in a region between adjacent mesh points. Although 
the electric field slope is assumed to be negligible, the electric field itself, 
in general, is allowed to vary from mesh point to mesh point. Moreover, under 
these conditions the remaining material properties also exhibit small changes 
between neighboring mesh points for low injection levels. These changes become 
large under high injection levels, in which case a more dense mesh point field 
may be required to maintain comparable changes. As a result, a near-exact 
closed-form solution is obtained at each mesh point for the photoexcited electron 
concentration, which is an exact solution of Eq. (9) and is represented by 



= B^.exp (w-yX) 


+ B^expCw^x) + Bjj (x) 


( 10 ) 


where the terms comprising Eq. (10) are given in Table 1. 

The subscripts, j, in Eqs. ( 10) — (30) repr isent the j^ mesh point in the p- 
region as indicated in Figure 1. Correspondingly, in Eq. (10), and B 2 

represent the constants of integration, and B^ in Eq. (11) represents the 

photoexcited electron concentration, produced by -ho ton absorption in the region 

of the j^k mesh point. In the context of the condition? 1 imposed on the continuity 
equations at each mesh point, and u 2 in Eq . (14) are constants. Photo- 
excited electron concentration, Eq. (10), is governed by he exponential tem.j 

and B 0 . which are all functions of position; in addition, B_. is also a function 
3j 3j 

of wavelengtn through Eqs. (12) and (13). Eq. (10) provides for continuous 

values of n^ e through the assignment of values to x, where the range of x is 

restricted to the region between ne (j-1) and j mesh points. In contrast to 
numerical integrate r methods, the continuum of values describing the behavior 
of n^ e is preserved in the closed-form scheme as was the original intent pro- 
posed by the Van Roosbroeck concept related to the use and the interpretation of 
the transport equations. 


It is clear that the relationship used in the closed-form aetnod, to 
represent the photoexcited electron concentration, Fq. (10), is an analytical 
solution that is demanded by, and has its support in, the transport equations. 
Moreover, this relationship may also incorporate a comprehensive set of phenomena 
submodels as dictated by the transport equations, imposed boundary conditions, 
and the representations c. the material properties coupled with the mesh point 
field. The phenomena submodels that influence Eq. '10) are represented in the 
relationships for the parameters, Eq. (11;-(3G), that describe their hehavicr in 
the region bounded by the (j-1) and j planes in Figure 1. Eqs. (IP _() describe 
the photon absorption generation rate and the redistribution of r 1 photoexcited 
electron-hole pairs, while Eq. (14) governs the electron effective diffusion 
length in the presence of an electric field. D’-ift and diffusion components of 
the electron current are represented in Eq. (lo), and the electric field used in 
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Table 1. List of relationships of terms appearing in t\>. solution of 
the electron concentration. Equation (10). 


3 j u ' 2 j ‘ w lj 


[G 1 ^ (x)exp(w^x) + G 2 .,(x)exp v 'J 2 jX)] 


(ID 


G lj (x) = 


i r x ^ v* * ) f*' 

— j e*p(- Ulj *’) j VV X P ( - J a j dX " 


)dAdx* ( 12 ) 


: 2j (x) = - d 1 : f exp( -“2j x ’) ( 

Jo J 


X (x’> 


a.N_.exp(- f a.dx’^dXdx 1 (13) 
303 ), 3 


“u 


2j 


buEi + // “"3 E pj \ | ~ - 

2 \j V l 2D nj / "J 


(14) 


dn . 

J . = qy E . + qD . 
nj nj ^3 p: n nj dx 


(15) 


E PJ " Vj + E pBNj + E plPj + E pPEj ’ 


(16) 


where 


b ^ J i 

? = PJ L- 

p£ij "Vpj 


(17) 


F = fh - £ 3 P.°j v 

pBNj S ; r pj pBNoj ’ 


(18) 


= _ kT 2_ dn iej 

pBNoj " q n ± dx 


(19) 
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Table 1 (continued) . 
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( 20 ) 


^21) 


( 22 ) 


( 23 ) 


( 24 ) 


( 25 ) 

( 26 ) 

( 27 ) 

( 28 ) 

( 29 ) 

( 30 ) 
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the drift term is given by Eq. (16). The electric field com nents are defined 
in Lqs. (17)-(27), arising under equilibrium and nonequilibrium conditions, and 
include effects produced by ohmic voltage drops in the quasi-neutral regions, 
bandgap narrowing, impurity concentration profile, and photoexcited electron 
concentration distribution. Eq. (29) describes the intrinsic concentration with 
respect to temperature and bandgap narrowing. 

In contrast, the relationship that is .^ed in rw methods, corresponding to 
Eq. (10), is typically linear and 's of the form given by 


n (x,, = n (x. . 
pe pe j-1 


n 

) + 


00 - V ( Vi > _ 

- x^ 1S (j-1)' 


?e‘ j 


(31) 


(j-1) 


where .. is the independent variable in the region bounded by the (j-1) and j 


planes. 


The photoexcited carrier concentrations, n . and n do not 

r pej pe^j+1) 


normally represent algebraic expressions, but represent discrete values and 
require initial estimates of the concentrations at their designated positions. 
The richness cf Eq. (10) and its associated parameters, defined in Eqs. (M)- 
(30), in representing device physics is clearly evident in contrast to the 
representation, in Eq. (31). This relationship, Eq. (31), is used in tt methods 
because it simplifies the matrix equation that requires solution, but it is 
artifical in its representation of the photoexcited electron cc acent ration 
because it has been constructed independent of the transport equations. 


The difference in tie results obtained by applying Eqs. (10) and (31) _o 
th same mesh point distribution, which defines the solar cell structure, is 
significant. Applying Eq. (10) results in the determination of the constants of 
integration, B^ and B 2j* ass lg ne£ ^ to the carrier analytical relationship at 

each mesh point. Substituting B.^ and B* . in Eq. (10) provides for an analy- 

tical relationship at each mesh point, describing the behavior of n .. f° r a 

continuum of values of x in the range x^ ^ and x^ . In contrast, the results 

of applying Eq. '?1) is a set of discrete values for n ^ at each mjsh. point. 

In the work reported in this paper, the general case is treated as it 
relates to injection level. The information of injection le /el is contained in 
the electric field and its components, Eos, (16)-(30), lifetime through the 
diffusion length, L^, and bounder’ conditions at the depletion edges bounding 

the p n junction, and in the carrier mobilities and dif fusivities . Application 
of p ->. (10) requires negligibly small changes in electric field between adjacent 
mesh points. Under these conditions, Eq. (10) is an excellent approximation to 
the exact solution . t w he assigned mesh point region. The measure for which 
Eq. (10) approaches tha txact solution is the sell-consistency obtained from the 
solution of the Poisson e , ration using the analytical relationships for n an d 

the numbe of iterations 1 'quired to obtain values of B, . and 

1.1 2j 

The fp-mesh point distribution in che p-type region o. the solar cell 

comprises (f 1) internal bor-^aries, and external boundaries at x., and x_ . At 
P J a 

the external boundaries, the usual p-n junction boundary condition on minority 
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carrier concentration applies. At each of the (f^-1) internal boundaries, the 

electron concentration and electron current, separately, are continuous and are 
represented by 


where 


n . (y .) * n , (0)k 4 

pej 7 oj pe(j+l) j 


k = e AEcj/kT 

j 


(32) 


(33) 


and 


J .(y .) = J ,..,*(0) 
nj 7 oj n(j+l) 


(34) 


where y is the separation between the (j-l) t ^ 1 and j t ^‘ mesh points. A surface 

recombination velocity boundary condition is imposed at the x,. boundary. This 

results in 2f^-boundary conditions, from which 2f^-equations are obtained that 

must be solved simultaneously to obtain the solution of the 2f -constants of 
integration. ^ 


Applying the boundary conditions using Eqs. (10) and (15) results in 2f - 
equations which are represented in matrix form by ^ 
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, and the 


The constants of integration are denoted B^, B ^ » ; » B 2 f 

relationships for the other factors contained in Eq. (35) a?e the P material and 
structure properties at corresponding mesh points. The Fletcher boundary 
condition [16] is used at x^, because it represents the boundary for the com- 


plete range of conditions that may exist: from equilibrium 

voltage V ^ = ^) through increasing positive values of V 

high injection levels where V 


Jph 


Jph 


^ V, 


bi' 


(i.e., the photo- 
up to and including 


In Eq (33), AE . represents the band edge discontinuity located at the j 
mesh point. If AE . is positive, the boundary condition describes minority 
carrier confinement (i.e., electrons confined to the junction region). 


The factors and u 2 j * ^q. (14), are the reciprocals of the effective 

electron dii fusion length in the p-type region, and describes the recombination 
of electrons in the presence of an electric field in the region bounded between 
the (j-1) and j mesh point. It describes the recombination related to those 
electrons entering this region across the planes defined by the (j-1) and j mesh 
points as well as those photoexcited carriers produced by photon absorption 

within that region. If the electric field, E . = 0, then w, . 

PJ ^ 

equal in magnitude to the reciprocal of the diffusion length, = /D^t^ 
Similarly, for degenerate material or for high injection levels, low values of 


- -<*>„ . and are 
2j. 


the electric field are obtained, and u) 


In both cases, electron drift 


Ij "2j ' 

toward the function occurs by means of diffusion rather than a combination of 
diffusion and field-assisted drift. Eq. (10) reduces to the more familiar form 
to represent photoexcited carriers. 


The constants of integration are obtained by solving Eq. (35). While the 
matrix in Eq. (35) may be inverted to obtain the solution, a recursion relation- 
ship exists between the constants of integration. As a result, there is a 
significant reduction in CPU execution time to obtain the values of these 
constants through the recursion relationships. For example, in certain iteration 
procedures, some of the terms in these relationships that depend on geometry and 
materials properties may not ch -nge and need not be calculated for every 
iteration. There are probably other benefits, which will be revealed as more 
experience is gained in using tnis type of modeling program. 


The recursion relationships for the constants of integration are given by 
the following equations: 



F (2f )- S n B 3£ ( >0f > - °f 6 nf < 6 f + V 
E. P ... £_ E E P 

T f (£ f + V - \f °f <S f + V 

P P P P P 


(36) 
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( 37 ) 
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B 2J 
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x.[e -6, .k .]-y .a. [6.-6, .k ] 
j j (i+I) nj nj ] j (j+1) nj 
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lj nj nj 2j 


(38) 

(39) 
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B 22 (6 2 £ 2 ) ^ F 01 6 2~ F 02 q l 9 nl ($ 2 k nl V 
T l (e r 6 2 k nl ) - Y nl a l (6 r 6 2 k nl ) 


(40) 


and 


B 11 = 9 n r Y nl B 21 


(41) 


There also exists recursion relationships for the 0 . ’s and the y .'s which are 

nj nj 

functions of material properties at their assigned mesh points. 


Having obtained the constants of integration and substituting them into the 
relationship representing the appropriate mesh point, there results f^- 

relationships which fully describe the behavior of n^ e in the p-region. In 

principle, all device properties may be recovered through the manipulation of 
n .s. Using Eqs. (10) and (15) and the corresponding relationships for holes 

in the n- region, the hole current J^, Eq. (28), is obtained. Eq. (28) represents 


a relationship relating photocurrent versus photovoltage. The photovoltaic J- 
Vjph curve ma y be obtained from which maximum power and efficiency is calcu- 
lated. Moreover, the effects on the J-Vj ^ curve of the phenomena submodels. 


and material and structure parameters may be investigated through B 


and B 


3j 


lj ’ 2j • 


SUMMARY 

A new computer modeling method is described and is applied to silicon solar 
cells. The method is similar to numerical integration (NI) methods in that 
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both require the use of a mesh point field. The set of transport equations 
is applied to each mesh point, and through a judicious selection of the mesh 
point distribution, an accurate closed-form analytical function is obtained at 
the assigned mesh point. Application of the boundary conditions, to an f-mesh 
point field, results in 2f-equations that require simultaneous solution. This 
solution is manifested through the determination of the 2f-constants of integra- 
tion, where each closed-form solution, representing a mesh point, contains two 
constants of integration. Solar cell transport solution is represented by the 
complete set of constants of integration obtained in the n- and p-regions. 
Substituting the constants of integration into the corresponding closed-form 
analytical function, representing an assigned mesh point, results in a set of 
analytical functions that is applicable only to its assigned region in the mesh 
point field. As a result, the complete set of closed-form functions describes 
the behavior of the photoexcited carrier concentration for a continuum of x- 
values in the n- and p-regions. The photoexcited carrier relationship is used 
to obtain the current-voltage relationship, from which the maximum power point 
and conversion efficiency are determined. Effects of temperature, solar concen- 
tration, submodel parameters, and structure parameters may be studied through 
changes at each mesh point in w-,., B^., f^., anc * B 3‘‘ 
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DISCUSSION 

(LAMORTE) 

LINDHOLM: I have a number of small questions. You put your attention on the 

dependent variables, which you treat as the hole and electron densities. 
I*m wondering how you enter the external boundary conditions since you 
are not feeding in hole and electron densities at the ends of the 
device. 

LAMORTE: At the junction side cf the end region, for example, meaning the 

emitter, the boundary condition becomes the Fletcher boundary condition 
for the photo-excited holes, because we are treating the general case of 
high injection level. At the surface region we are calculating the re- 
combination velocity, by means of a Ax divided by the lifetime at the 
center of the first mesh point. And we have gotten some interesting 
results there, in that it appears that even for calculating the surface 
recombination velocity in that manner, it may not be consistent with the 
slope of the photo-excited holes at the surface. If you don't recal- 
culate the surface recombi nation velocity, the solution will oscillate. 
We're just putting a fix in that, and the fix is that you want to update 
the surface recombination velocity at the surface by the exact relation- 
ship. Meaning the product of the surface recombination velocity times 
the photo-excited hole concentration at the surface, and that is equal to 
the surface recombination current, which includes a drift field component 
and a diffusion component. So that becomes the left-hand boundary. 

LINDHOLM: I have a related question. If you are going to do a non- illuminated 

analysis where you apply a voltage, then bow do you get into the external 
boundary conditions? I'm not now concerned with the edge of what you call 
the depletion region but rather the contacts of the device. 

LAMORTE: I don't know whether it will work for that, but we c-n probably make 

it work. We have not given that any consideration. 

LINDHOLM: Would you integrate .he electric field? Would that be the way of 

doing it? Getting integrated, the electric field through the material, 
and setting that equal to zero, and then you would have to iterate, I 
suppose . 

LAMORTE: Yes, we have that in the model because under high injection level 

you want to determine what the voltage drop is across the quasi- neutral 
region. 

LINDHOLM: When you say the Fletcher boundary condition, you mean the Fletcher 

Masawa boundary condition as modified by Houser? 

LAMORTE: No, the Fletcher and Masawa conditions are separate. They account 

for the same thing, but they require different information. The Fletcher 
boundary condition applies to the edges of the depletion region and thats 
what wo are using. The Masawa uses the right-hand edge of the depletion 
region and the left-hand contact. 
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LINDHOLM: How do you define the edge of the depletion region in your numerical 

analysis? 

LAMORTE: We calculate the depletion region for either a graded case or the 

abrupt case. 

LINDHOLM: You calculate it how? Depletion approximation? 

LAMORTE: Yes. 

LINDHOLM: That won't work, I'm afraid, because you probably need to make a 

correction because of the free electrons and holes in that region. We 
could discuss that privately. That is probably an updating of the physics 
that goes into the model. 

LAMORTE: Well, that is updated. 

LINDHOLM: Oh, I see. You don't use the depletion approximation to calculate 

the depletion thickness? 

LAMORTE: We use it to get it started. 

LINDHOLM: Oh, I see, then you update. Okay. Then how do you define the edge 

of the depletion region once you get it updated? Somehow you have to use 
some criterion to define the edge of the depletion region. 

LAMORTE: Then we use the depletion approximation. 

LINDHOLM: Okay, I don't understand the answer to that, but maybe we could 

discuss that privately, unless you want to elaborate on it now. What I 
thought you said was, as your starting point to get the edge of what you 
call the depletion region, you would use the depletion approximation and 
then you would iterate up in, including the electron hole densities, in 
order to redefine the depletion approximation. Then I asked you what the 
criterion was after you did the updating. 

QUESTION: Maybe that's something the two of you could discuss. 

LINDHOLM: I have one last question. This is a very interesting idea to me. I 
have been sounding very negative; I'm sorry to sound negative. I just 
was looking at some small points. I don't do numerical analysis, but I'm 
somewhat familiar with what Mike did, and Gummel, and all these other 
people, and the only place I have ever seen this done previously, similar 
to this, is in a book by some Russians. I wonder if this is the first 
time this method has been used? 

LAMORTE: I have not seen it anywhere else, and I hive spoken to about a dozen 
people who have done computer modeling in *>th< r "re as , other than semi- 
conductors, as well, and they claim that they nave not seen it. And I 
have not seen it elsewhere. 

LINDHOLM: This Russian book is not quite numerical analysis, so I guess I 

haven't seen it either, but that's as close as I have come to seeing it. 
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It * 8 very interesting. 

LAMORTK : One of the things that this does for you, and what we showed in the 

avalanche photo diode -- because you're reflecting a lot of the physics 
in the closed form solution -- when you look at the convergence of the 
solution as a function of mesh points, you find that you do not need as 
many mesh points, for example, in the avalanche photo diode. We didn't 
need to go above 15 to 20 mesh points when we had a half-micrometer 
depletion width, and we went all the way up to 80 mesh points, and at 15 
to 20 we were within IX of the asymptotic solution with 80 mesh points. 

I'm fully convinced that since you have a closed-form solution, which is 
a good approximation, the physical system, that it's forgiving in term.', 
of using a lesser number of mesh points. And so, therefore, that combined 
with the recursion formulas, the CPU time is reduced greatly. 

LINDHOLM: I was wondering if you could describe the output that you got from 

your model and do you have it set up to give you graphs of carrier dis- 
tributions and so forth? 

LAMORTE: No, we haven't done that yet. We almost have the model working, 
meaning that with one of your cells we got something like 38 milli- 
amperes, and it's about 20X too high, so we have gotten it that far and 
we are still trying to debug and determine where that is. And it may be 
just a simple matter that the lifetime we are using is too high. 

LIMDHOLM: So the code doesn't give you plots of carrier distributions and 

things that help one see what's going on in the physics of the device. 

LAHORTE: We get discrete points, for example, of carrier concentration, okay, 

at the mesh points, and we get it at two points actually, we get it at 
each of the mesh points and in between. 6ut if you wanted to, by taking 
the relationship that applies to that particular region of the cell in 
combination with the appropriate constants of integration, you could plot 
the entire thing on the continuum if you wished. 


369 



